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Abstract. We apply the moment analysis technique to analyze large scale simulations of the Zhang sandpile 
model. We find that this model shows different scaling behavior depending on the update mechanism used. 
With the standard parallel updating, the Zhang model violates the finite-size scaling hypothesis, and it 
also appears to be incompatible with the more general multifractal scaling form. This makes impossible 
its affiliation to any one of the known universality classes of sandpile models. With sequential updating, it 
shows scaling for the size and area distribution. The introduction of stochasticity into the toppling rules 
of the parallel Zhang model leads to a scaling behavior compatible with the Manna universality class. 

PACS. 05.70.Ln Nonequilibrium and irreversible thermodynamics - 05.65.+b Self-organized systems 



The identification of universality classes is one of the 
most important and still open problems in the field of 
self-organized criticality (SOC) [lj. In spite of the great 
relevance of this issue, however, it has not been possible 
until very recently to clearly discern differences in the crit- 
ical behavior of the various SOC models proposed so far. 
The situation seems now to have been settled in the case 
of the Bak-Tang-Wiesenfeld (BTW) § and the Abelian 
Manna models that represent, respectively, the proto- 
typical examples of deterministic and stochastic sandpile 
automata. In this particular case, recent large scale sim- 
ulations ||[||^,|^| clearly indicate that Manna and BTW 
models belong to different universality classes. 

The universality class asset remains still uncertain for 
the Zhang model ||, which is the archetype of all sand- 
pile automata with continuous variables. This model has 
deterministic dynamics like the BTW model, and in con- 
trast with most other cases, it is non-Abelian. This means 
that the stable configurations obtained at the end of an 
avalanche depend on the order in which the active sites 
are updated. In spite of this essential peculiarity, earlier 
simulations of the model || placed it in the same uni- 
versality class as the BTW model, of which it was sup- 
posed to be the continuous counterpart. This conclusion 
was confirmed by the large scale simulations performed by 
Liibeck in d = 2 [[To) . On the other hand, Milshtein et al. 
||, analyzing different magnitudes than Liibeck, arrived 
at the opposite result in d = 2, namely, they observed no- 
ticeable differences in the critical exponents. Finally, the 
simulations of Giacometti and Diaz-Guilera [jll| provided 
evidence that, even though the exponents of both models 
are similar in d — 2, they do not coincide in d = 3. 

Recently, it has been shown JlJ,|l3| that the determin- 
istic nature of the BTW model renders its scaling incom- 



patible with the standard finite-size scaling (FSS) hypoth- 
esis, and induces moreover peculiar non-ergodic effects ||. 
Thus, it would not be surprising to find similar anomalies 
in the Zhang model because of its deterministic nature. In 
this paper we characterize the Zhang universality class by 
applyi ng the recently proposed moment analysis technique 
[12lO]. In the following we will show that the scaling of 
the Zhang model depends very strongly on the updating 
mechanism implemented in the simulations, either paral- 
lel or sequential. The Zhang model with parallel updating, 
as it has been customarily defined in the literature, dis- 
plays a complex scaling behavior that is not compatible 
neither with the standard finite-size scaling (FSS) hypoth- 
esis nor with the multifractal picture Q proposed for the 
avalanche distribution of the BTW model Jl^|. On the 
other hand, the Zhang model with sequential updating 
shows well defined size and area exponents. The origin 
of the complex behavior of the parallel Zhang model can 
be ascribed to the deterministic microscopic dynamical 
rules of the model. In order to prove this conjecture, we 
study a variation of the Zhang model, the stochastic par- 
allel Zhang model, which exhibits a standard FSS behav- 
ior, compatible with the universality class of the Manna 
model. 

We consider the definition of the original Zhang model 
[^). On each site of a d dimensional hypercubic lattice of 
size L we assign a continuous variable Ei, called "energy" . 
At each time step, an amount of energy i5 is added to a 
randomly selected site j (Ej — > Ej+S) . The quantity 5 is a 
random variable uniformly distributed in [0, (5 max ] |lq ] . In 
our simulations we consider the fixed value <5 max = 0.25. 
When a site acquires an energy larger than or equal to 
lj (Ek > 1), h becomes active and topples. An active 
site k relaxes losing all its energy, which is equally dis- 
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tributed among its nearest neighbors: Ek — > 0, and E^ — > 
Ey + Ek/2d. Here the index k' runs over the set of all 
nearest neighbors of the site k. The transported energy 
can activate the nearest neighbor sites and thus create an 
avalanche. Energy can be lost only at the boundary of the 
system (open boundary conditions). The avalanche stops 
when all sites in the lattice are subcritical (Ei < 1). Given 
these dynamical rules, it is easy to see that the Abelian 
nature of the model depends on the type of updating im- 
plemented. With parallel updating — parallel Zhang (P-Z) 
model — at each time step t in the evolution, all active 
sites arc toppled simultaneously, and time is incremented 
t — > t + 1. Since all the energy of an active site is trans- 
ferred to its nearest neighbors, we notice that in a bipar- 
tite lattice (such as the hypercubic lattice here considered) 
all active sites at a give time step t reside onto the same 
sublattice, and that activity alternates between the dual 
sublattices in consecutive time steps. Again, since all the 
energy is transferred in a toppling, the state of the active 
sites in a sublattice at time t is independent of the or- 
der in which the active sites in the dual sublattice were 
updated at time t — 1. On the other hand, with sequen- 
tial updating — sequential Zhang (S-Z) model — at each 
time step t an active site is randomly chosen among all 
the N a (t) active sites present at that time. The chosen 
site is toppled and time is incremented t — > t + l/N a (t). 
In this case, activity is not restricted to alternate sublat- 
tices, but spreads all over the system. Depending on the 
order in which the intermediate list of over-critical sites is 
updated, any active site with at least one active nearest 
neighbor can end up with an energy Ei — (if its nearest 
neighbors are updated before it) or with energy Ei > 
(if it is updated before its nearest neighbors). In this case 
the model fully exploits its non- abelian characters. 

We have analyzed both parallel and sequential Zhang 
models by determining their critical exponents jlj] . In the 
limit of an infinite slow driving, i.e. the energy addition 
is interrupted during the avalanche evolution, defining a 
complete time scale separation the system reaches 
a critical stationary state with avalanches of activity dis- 
tributed according to power laws. If we define the proba- 
bility distributions P(x) of occurrence of an avalanche of 
a given size s, time t, and area a, the FSS hypothesis fll7| , 
usually assumed in SOC systems, states that 



P{x,L) 



(i) 



with x — s,t, or a, respectively. If the FSS ansatz is valid, 
then the critical exponents (5 X and t x completely deter- 
mine the universality class of the model under scrutiny 

Previous numerical works on the Zhang model P,[10l 
[TT[ | have most often proceeded measuring the exponents as 
the slope in a log-log plot of the density P(x, L) as a func- 
tion of the magnitude x. Even though with this procedure 
one can determine the exponents within a 10% accuracy, 
its performance is affected by the existence of the upper 
and lower cutoffs, which render difficult its application. In 
this respect, it is better to use analysis techniques that 
contain explicitly the system-size dependency. 



The moment analysis technique was introduced by De 
Menech et al. in the context of the two dimensional BTW 
model [Ell, and its validity has been extensively checked 
for both Abelian and stochastic models H,fjJ. In this me- 
thod, the q-t\i moment of a probability distribution on 
a lattice of size L is defined by (x q Yr — J x q P(x, L)dx. 
Assuming the FSS hypothesis, Eq. (pi), the g-th moment 
has the size dependence: 



* T{y)dy ~ L^ q+1 ~ T '\ 



(2) 

where we have introduced the transformation y = x/L^ x . 
In general, one has (x q ) L ~ L^W, where the exponent 
a x (q) can be obtained as the slope of the log- log plot of 
(x q ) L as a function of L. If the FSS hypothesis is indeed 
correct, we expect u x {q) ~ fi x {q + 1 — t x ), and therefore 
one can compute (3 X — da x (q)/dq. For very small values 
of q this is not correct, since the integral in (||) is dom- 
inated by its lower cut-off. Once computed the exponent 
(3 X , the corresponding t x is obtained using the scaling re- 
lation <7 X (1) = f3 x (2 - T x ). 

In order to apply the moment analysis technique, we 
have performed simulations of the P-Z and S-Z models in 
d = 2, for sizes ranging from L = 128 to L = 1024. Statis- 
tical distributions are obtained by averaging over 5 x 10 6 
nonzero avalanches. As a consistency check of our algo- 
rithm, we have estimated the average energy of the sys- 
tem in the stationary state, defined by E — (J^i Ei) /L 2 , 
where the brackets denote an average over at least 10 6 
stable configurations. Extrapolating to an infinite system 
size, we obtain a value E = 0.630 ± 0.005 for the P-Z 
model and E = 0.626 ± 0.005 for the S-Z model, both in 
good agreement with previous estimates (E ~ 0.62 — 0.63) 
10 ■ 

We have computed the moments cr x (q) for our data 
from both the P-Z and S-Z models. In the presence of 
FSS, as we have argued above, one should expect to ob- 
serve da x {q) / dq = [3 X = const. Simulations on the Manna 
model, Ref. 0, show that this is indeed the case, with 
a slope for the moments a x {q) that reaches a saturation 
value for relatively small values of q JTsf ] . 

In Fig. [j] we have plotted da x (q) / dq for the P-Z model 
as a function of q, for q < 2.5. We observe that all the 
moments present a noticeable curvature for all q, and do 
not seem to reach a constant slope for the values of q 
considered. For values of q larger than 4, the slope of the 
moment functions seems to finally achieve a saturation 
value. As an estimate of this trend, in Table jl we report 
the value of the local slope at different q's, [13|, defined 
by Aa x (q) — cr x (q + 1) — cr x (q). Assuming that the cur- 
vature of the moments at small q is merely a crossover 
effect and that the real exponents f3 x are given by the 
saturation plateaus at large q, we are led to the values 
(3 S = 3.39, t = 1.74, and f3 a = 2.25. In particular, the 
value of P a is completely unphysical: The maximum area 
of an avalanche cannot grow faster than L d in d dimen- 
sions, and thus f3 a must be smaller than 2. This tells us 
that the FSS form used in Eq. (||) is not adequate in the 
case of the P-Z model, and leads to spurious results. In 



Romualdo Pastor-Satorras, Alessandro Vespignani: Anomalous scaling in the Zhang model 



3 




Fig. 1. Derivative of the exponents a x (q) for the parallel Zhang 
model. The monotonous increase of the exponent indicates the 
lack of scaling in this model. 




Fig. 2. Data collapse analysis of the integrated distribution 
of sizes, for the parallel Zhang model. System sizes are L — 
128, 256, 512, and 1024. Exponents are fj s = 3.39 and r s = 1.42. 



what respects to the size and time distributions, we can 
check the validity of this result by means of a data collapse 
technique: If the FSS ansatz Eq. (0) is correct, then by 
rescaling x — > x/L^ x and P{x) — > L^ xTx P(a;), we should 
obtain distributions that collapse onto the same universal 
curve for different values of L. In Fig. ^ we plot the data 
collapse of the size distribution P(s), with the exponents 
fi s = 3.39 and r s = 1.42,. The really poor collapse of the 
curves testifies that also the avalanche size distribution 
does not fulfill the FSS in the P-Z model. We observe a 
similar lack of collapse for the area distribution. 

From Figs. and |, and Table we conclude that the 
Zhang model with parallel updating violates the FSS hy- 
pothesis, Eq. (||). Noticeably, this lack of FSS is observed 
also in numerical simulations of the directed version of 

g 1 2 3 4 5 

Aa s (q) 3l7 330 337 339 OfT 

Aat(q) 1.29 1.61 1.70 1.73 1.74 

Aa a (q) 2.09 2.19 2.23 2.25 2.26 

Table 1. Local slope of the moment exponents a x (q) in the 
parallel Zhang model, for different values of q. 
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Fig. 3. Data collapse analysis of the integrated distribution 
of sizes for the sequential Zhang model. System sizes are L — 
128, 256, 512, and 1024. Exponents are r s = 1.29 and f3 s = 2.78. 

the Zhang model |^(|. We have also tried to fit our data 
to the more general multifractal scaling form proposed by 
Kadanoff et al. [H, and applied to the BTW model in 
Ref. In this form of multifractal analysis, one tries 
to collapse the data to the form log(P(x, L))/ log(aL) as 
a function of log(x)/log(aL), for a suitably chosen con- 
stant a. We have checked that this sort of scaling is also 
not compatible with our data of the original Zhang model. 
In particular, we have not succeeded in finding a constant 
a for which the scaling is correct. 

The moment analysis of the S-Z model yields good re- 
sults for the size and area distributions, with derivatives 
of the moments cr x (q) that reach a saturation plateau for 
small values of q. The values obtained are r s = 1.29(2) and 
P s = 2.78(2) for the size exponents, and in Fig. || we plot 
the data collapse analysis for the size distribution. The 
perfect collapse of this figure confirms the good scaling of 
this model. The same is obtained with the area distribu- 
tions with exponents r Q = 1.43(2) and [3 S = 1.94(2). The 
time distribution, however, shows a lack of scaling, with 
a not clear trend for the derivative fit = da t (q)/dq as a 
function of q. This feature could have several origins in- 
cluded the possibility that in this case we do not have yet 
reached the scaling regime because of the different time 
updating, that gives a very small range of time scales. 

The complex quality of scaling in the parallel Zhang 
model can be attributed to the deterministic nature of 
its dynamical rules, which is somehow smoothed by the 
stochastic updating in the sequential model. In order to 
check this conjecture, we propose a variant of the P-Z 
model, the stochastic parallel Zhang model (SP-Z), in 
which energy is stochastically redistributed. The model 
is defined by the following modifications of the relaxation 
rules. In the SP-Z model, an active site k loses all its en- 
ergy, Ej~ — > 0, which is randomly redistributed among its 
nearest neighbors. In the practical implementation of this 
rule, we draw four random numbers < < 1, 
with ^2 k , £k' = 1 and update the nearest neighbors k' 
by Ek' — * Ek' + £k'Ek plj] . In this model, the stochastic 
energy input is a random variable uniformly distributed 
in [0, <5 max ], with 5 max = 0.25. Sites still have a contin- 
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Fig. 4. Data collapse analysis of the integrated distribution of 
sizes for the stochastic parallel Zhang model. System sizes are 
L = 128, 256, 512, and 1024. Exponents from Table §. 



uous spectrum of energy, but the new dynamical rules 
are stochastic. If the assumption that stochasticity yields 
a standard scaling behavior, then this model should be 
regarded as the continuous counterpart of the original 
Manna model § . 

We have performed numerical simulations of the SP- 
Z model in system sizes ranging from L = 128 to L = 
1024, averaging the probability distributions over 5 x 10 6 
nonzero avalanches. We observe that slope of cr x (q) reaches 
a saturation value for very small values of q, for sizes, 
areas, and also times. In Table || we report the values 
obtained for the exponents j3 x and t x . As expected, the 
exponents are in perfect agreement with the values found 
in the Manna model. This fact confirms the presence of 
a unique universality class for all stochastic models. As a 
final check, we show in Fig. |^ the data collapse analysis 
for the distribution of sizes. The perfect collapse of these 
plots should be compared with the poor result obtained 
in the original Zhang model, shown in Fig. ^. 

In summary, applying the moment analysis technique, 
we have shown that the scaling of the Zhang model de- 
pends on the updating rules implemented in the simula- 
tions. The Zhang model with parallel update violates the 
FSS hypothesis Eq. ([[]) for the avalanche distributions of 
sizes, times, and areas. The anomalous scaling is stronger 
than in the BTW model, since data cannot be fitted even 
to the more general multifractal scaling form. This, con- 
trary to previous claims |^, |Io| , pT[ , makes impossible to 
assign any precise universality class to this model. The se- 
quential updating introduces a small amount of stochastic- 
ity in the Zhang model that, in this case, shows FSS for the 
size and area distributions. In spite of this property, the 



lack of scaling for the time distribution does not allow to 
place the sequential Zhang model in a definite universality 
class. The anomalous scaling of the original Zhang model 
can be therefore ascribed to the deterministic nature of the 
dynamical rules defining the model. A stochastic versions 
of the model show a standard FSS behavior, compatible 
with the universality class of the Manna model. 

This work has been supported by the European Net- 
work under Contract No. ERBFMRXCT980183. We thank 
A. Stella, A. Vazquez, and S. Zapperi for helpful com- 
ments and discussions. 
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